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ABSTRACT 

Context. Since the end of 2013 a new generation of instruments optimized to image young giant planets around nearby stars di¬ 
rectly is becoming available on 8-m class telescopes, both at Very Large Telescope and Gemini in the southern hemisphere. Beyond 
the achievement of high contrast and the discovery capability, these instruments are designed to obtain photometric and spectral 
information to characterize the atmospheres of these planets. 

Aims. We aim to interpret future photometric and spectral measurements from these instruments, in terms of physical parameters of 
the planets, with an atmospheric model using a minimal number of assumptions and parameters. 

Methods. We developed Exoplanet Radiative-convective Equilibrium Model (Exo-REM) to analyze the photometric and spectro¬ 
scopic data of directly imaged planets. The input parameters are a planet’s surface gravity (g), effective temperature (T^jf), and 
elemental composition. The model predicts the equilibrium temperature profile and mixing ratio profiles of the most important gases. 
Opacity sources include the H2-He collision-induced absorption and molecular lines from eight compounds (including CH4 updated 
with the Exomol line list). Absorption by iron and silicate cloud particles is added above the expected condensation levels with a fixed 
scale height and a given optical depth at some reference wavelength. Scattering was not included at this stage. 

Results. We applied Exo-REM to photometric and spectral observations of the planet /? Pictoris b obtained in a series of near-IR filters. 
We derived Aff = 1550 ± 150 K, log(g) = 3.5 ± 1, and radius R = 1.76 ± 0.24 Rj^p (2-cr error bars from photometric measurements). 
These values are comparable to those found in the literature, although with more conservative error bars, consistent with the model 
accuracy. We were able to reproduce, within error bars, the J- and H-band spectra of [} Pictoris b. We finally investigated the precision 
to which the above parameters can be constrained from SPHERE measurements using different sets of near-IR filters as well as 
low-resolution spectroscopy. 

Key words, radiative transfer, planets and satellites: atmospheres, planets and satellites: gaseous planets, stars: individual (/? Pictoris) 


1. Introduction 

Following the detection of 51 Peg b by Mayor & Queloz (1995) 
using velocimetry, almost 2000 exoplanets' are known as of 
today (March 10, 2015) and many more candidates are awaiting 
confirmation (Rowe et al. 2015). Among them, only a few were 
detected with direct imaging. 

The first image of a planetary mass object orbiting a star, 
2M1207 b, was obtained by Chauvin et al. (2004) with NaCo 
at the Very Large Telescope (VLT) and this result has inspired 
several other discoveries in the last decade (Marois et al. 2008; 
Lagrange et al. 2009; Rameau et al. 2013). In the first case, the 
mass ratio was highly favorable as the central star is a brown 
dwarf (BD). The detection was enabled by the use of an adaptive 
optics (AO) system in the L’ band and no specific device to 
attenuate the star, such as a coronagraph, was needed. Later, 
larger mass ratios became feasible with the improvement of 
high contrast imaging techniques. For now, the planet with the 


largest mass ratio with respect to its host star, and for which we 
have an image, is HD95086 b with a mass of 5 + 2 Mjup around 
a star of 1.6 Msun (Rameau et al. 2013). Conveniently, direct 
imaging also allows us to collect spectroscopic data if one is 
able to attenuate the starlight at the location of the planet. Janson 
et al. (2010) presented the first spatially resolved spectra of 
HR8799 c still with NaCo. The spectrum covered the 3.88-4.10 
pm interval but with a low signal-to-noise ratio. The OSIRIS 
instrument at Keck II allowed Bowler et al. (2010) and Barman 
et al. (2011) to obtain spectra of the exoplanet HR8799 b in the 
K band, sensitive to methane opacity, as well as a spectrum in 
the H band (Barman et al. 2011). Later, the same instrument 
provided spectra of the same planet Barman et al. (2015) and 
of HR8799 c (Konopacky et al. 2013) in the K band at higher 
resolution (R ~ 4000). Near-infrared, low-resolution spectra of 
the planet fi Pictoris b (Lagrange et al. 2010) were obtained with 
GPI (Macintosh et al. 2014) in the J and H bands (Bonnefoy 
et al. 2014; Chilcote et al. 2015), an instrument tailored for the 
search of young giant planets. 


' exoplanet.eu 
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Fig. 1. Spectra calculated for a cloudy (blue) and cloud-free (red) atmo¬ 
spheres at a distance of 10 pc, for a radius of 1 log(.g) = 4 and Feff 
= 1400 K . Examples of BT-Settl (cyan) and Drift-PHOENIX (gray) 
models with the same parameters are also shown for comparison. 

At the present time, instruments for direct imaging combine 
AO, coronography, and differential imaging to detect faint 
planets. The Spectro-Polarimetric High-contrast Exoplanet 
Research, SPHERE (Beuzit et al. 2008), installed at the VET, is 
designed to perform high contrast imaging for detecting young 
giant planets and for characterizing their atmospheres. SPHERE 
provides broad- and narrowband photometry and spectroscopy 
in the near-infrared (NIR) range with the InfraRed Dual-band 
Imager and Spectrograph (IRDIS, Langlois et al. 2010) and 
the Integral Eield Spectrograph (lES, Claudi et al. 2008), and 
photometry and polarimetry in the visible range with the Zurich 
Imaging Polarimeter (ZIMPOL, Schmid et al. 2010). 

To directly detect the light from a planet around a star other 
than the sun, the following conditions have to be met; 

- the star-to-planet angular separation must be larger than the 
angular resolution offered by an 8 -m telescope in the NIR 
(25-50 mas for SPHERE). This restrains the sample of tar¬ 
gets to less than 100 pc as well as the minimal physical sep¬ 
aration to >1 AU. 

- the star-to-planet brightness ratio must be smaller than the 
achievable instrumental contrast, which is typically 10 ^- 10 ^ 
at less than 1". Only giant planets can be warm enough at 
young ages to produce a detectable emission (Burrows et al. 
1995; Chabrier et al. 2000). These young extrasolar giant 
planets (YEGP) are ~ 10-100 millions years old. 

One theoretical challenge is to understand planetary forma¬ 
tion mechanisms. Brown dwarf and young extrasolar giant planet 
are two types of objects almost impossible to differentiate, which 
are formed in different ways. Proposed mechanisms are gravita¬ 
tional instabilities (Boss 2001), such as stars, without a suffi¬ 
cient mass to start to burn hydrogen or core accretion (Lin & Ida 
1997). A precise determination of the luminosity, mass, and age 
can supply information about the initial entropy of the planet and 
allows us to identify the formation mechanism (Bonnefoy et al. 
2014; Marleau & Gumming 2014). 

In parallel, atmospheric models for objects with mass and 
temperature lower than an M dwarf were developed since the 
end of the 1990’s. The basic idea is to include some chemistry 
and other physical processes in a H-He atmosphere to account 
for the range of pressure temperature expected in such low-mass 
objects. The models differ in particular by their treatment of 
dust opacity. This is a crucial component of the models since 



Fig. 2. Spectra calculated for a cloud-free atmosphere at a distance of 
10 pc, for a radius of 1 Rj„p, log(.g) = 4.0 and two different effective tem¬ 
peratures (Teff): 700 K (blue) and 2000 K (red). Molecular absorptions 
other than H 2 O are indicated. Examples of BT-cond models (cyan and 
gray) are also shown for comparison. 

it was noticed very soon that these atmospheres must contain 
clouds below T^-fi ^ 2600 K to account for the spectroscopic 
observations (Tsuji et al. 1996). 

The model of Tsuji and collaborators is a direct adaptation of 
their M dwarf models to BD. Dust grains in the atmosphere are 
treated through a parametrized model (Tsuji 2002), with three 
cases. In Case B, the part of the atmosphere where the temper¬ 
ature is lower than the condensation temperature is full of dust. 
In an opposite case (Case C), dust forms but is immediately re¬ 
moved by precipitation so that it does not contribute to the opac¬ 
ity. The third case, called the Unified Cloud Model, is intermedi¬ 
ate between these two extremes with dust present only between 
the condensation level and a level at a slightly lower temperature 
Tc. 


Table 1. Atom and molecular opacity sources 


Opacities 

Intensity cutoff 
(cm molecule"^) 

References 

H 2 O 

CO 

10^^' at 2500 K 

10-27 at 3000 K 

HITEMP line list 

(Rothman et al. 2010) 

CH 4 

10-27 at 1500 K 

Yurchenko & Tennyson (2014) 

NH 3 

5 X 10-27 at 4000 K 

Yurchenko et al. (2011) 

TiO, VO 

10-22 at 4000 K 

Plez (1998) (with update) 

Na, K 


Kramida et al. (2014) 
Burrows & Volobuyev (2003) 

H 2 -H 2 , 

H 2 -He 


Boiysow et al. (1988, 1989) 
Borysow & Frommhold (1989) 
Borysow et al. (2001) 
Borysow (2002) 


Marley et al. (2000) proposed a model adapted from solar 
system planets and using the Ackerman & Marley (2001) cloud 
model, which is parametrized with a factor representing the 
sedimentation efficiency. 

The Lyon’s group atmospheric models, DUSTY and COND 
(Allard et al. 2001) are two cases similar to Tsuji’s B and 
C cases. BT-Settl (Allard et al. 2003) is a more complex 
model, which compares condensation, sedimentation, and mix¬ 
ing timescales of dust to define cloud parameters. 
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Finally, Drift-PHOENIX (Helling et al. 2008) consid¬ 
ers microphysical processes (nucleation, condensation, particle 
growth, sedimentation, and evaporation) to calculate composi¬ 
tion, number density, and size distribution of dust particles as a 
function of atmospheric pressure level. 

All these models were used to constrain the main parameters, 
such as the effective temperature, surface gravity, atmospheric 
compounds, or radius of the detected planets (Neuhauser et al. 
2005; Marley et al. 2012; Bonnefoy et al. 2013; Galicher et al. 
2014; Bonnefoy et al. 2014) and to predict abilities of new 
instruments (Boccaletti et al. 2005; Hanot et al. 2010; Vigan 
et al. 2010). Most of these models are developed for BD and 
applied to exoplanets. Although, the radius of exoplanets can 
be similar to that of BD, they have lower masses. Therefore, 
the range of surface gravity considered in BD atmospheric 
models (log(g)>3.5) does not necessarily cover the entire range 
expectable for YEGP (log(g)>2). 

Direct imaging of YEGP is characterized by low flux, low 
signal to noise and low spectral resolution. In that respect, the 
models that are used to interpret these images should be repre¬ 
sentative of the level of data quality. Eor that purpose we specif¬ 
ically developed a model to analyze direct imaging of YEGP 
for instruments like SPHERE. It is a radiative-convective equi¬ 
librium model, assuming thermochemical equilibrium for self - 
luminous planets in which stellar heating is neglected. It allows 
us to explore low surface gravity, i.e. low-mass YEGP. 

The radiative-convective equilibrium model is described in 
Section 2. In Section 3, we apply Exo-REM to the well-known 
planet fi Pictoris b, derive physical parameters from existing 
measurements, and compare our results to previously published 
investigations. In Section 4, we analyze the uncertainties in the 
derived physical parameters as a function of photometric errors 
in the context of SPHERE observations. The conclusion is drawn 
in Section 5. 


2. Model description 

2.1. Radiative-convective equiiibrium modei 

2.1.1. Numerical method 

We solve for radiative-convective equilibrium, assuming that the 
net flux (radiative H- convective) is conservative, and neglecting 
the stellar flux impinging on the planet. The net fluxes are cal¬ 
culated between 20 cm ' and 16000 cm ' over 20 -cm ' inter¬ 
vals using the radiative transfer equation with no scattering and 
a k-correlated distribution method to represent molecular opac¬ 
ity. The atmospheric grid consists of 64 logarithmically equally 
spaced pressure levels. The system of equations fixing the con¬ 
stancy of the net flux over the atmospheric grid is solved itera¬ 
tively through a constrained linear inversion method. Details are 
given in Appendix A. 

2.1.2. Spectroscopic data 

As mentioned above, the spectral flux was calculated over 20 
cm ' intervals using a k-correlated distribution method. Eor each 
molecule and each interval, we calculated a set of Uk = 16 k- 
coefficients (l-l, Uk), 8 for the interval [0:0.95] of the normal¬ 
ized frequency g*, and 8 for the interval [0.95:1.00]. The val¬ 
ues of g* and associated weights m\ are those of the 8 -point 
Gaussian-Legendre quadrature for each of the two g* intervals. 
The k-coefhcients were calculated for a set of 15 pressures be¬ 


tween 100 bar and 0.01 mbar (2 values per decade) and, for each 
pressure, a set of 6 temperatures, increasing with pressure to en¬ 
compass model temperature profiles encountered in the literature 
for exoplanets with 500 K < Teff < 2000 K. Absorptivity spec¬ 
tra for a given pressure and temperature were calculated using 
a line-by-line radiative transfer program with a frequency step 
equal to the Doppler half-width of the lines. 

We considered the eight most important molecules and atoms 
in terms of opacity for relatively cool exoplanets (500 K < < 

2000 K): H 2 O, CO, CH 4 , NH 3 , TiO, VO, Na, and K. The origin 
of the line lists and the intensity cutoff used to calculate the ab¬ 
sorptivity spectra are given in Table 1. In the previous version of 
the model used by Baudino et al. (2013, 2014a,b), Galicher et al. 
(2014), and Bonnefoy et al. (2014), the methane line list origi¬ 
nated from Albert et al. (2009), Boudon et al. (2006), Daumont 
et al. (2013), and Campargue et al. (2012) for CH 4 , and from 
Nikitin et al. (2002, 2006, 2013) for CH 3 D. Our new methane 
line list now comes from the Exomol database (Yurchenko & 
Tennyson 2014). 

Table 2. Compounds considered in thermochemical equilibrium calcu¬ 
lations 


Compounds of interest' 

Species included in chemical 
equilibrium calculations' 

H 2 O, CO, CH 4 

H 2 O, H 2 O*, CO, CH 4 

NH 3 

NH 3 , N 2 , NH 4 SH*, H 2 S 

Na,K 

Na, Na 2 S*, H 2 S, HCl, NaCl, 

K, KCl, KCl*, NH 3 , NH 4 CI* 

TiO, VO 

Ti, TiO, Ti 02 , V, VO, VO 2 , 

Ca, CaTi 03 *, VO*, H 2 O 

Mg2Si04*, MgSi 03 *, Si02* 

Mg, SiO, H 2 O, Mg 2 Si 04 *, 
MgSi03*, Si02* 

Ee* 

Ee, Ee* 


': Species marked with asterisks are condensates 


Eor all species except alkali, we calculated line absorption up 
to 120 cm ' from line center using a Voigt profile multiplied by a 
X factor to account for sub-Lorentzian far wings. Foix, we used 
the profile derived by Hartmann et al. (2002) for H 2 -broadened 
lines of methane. The far wing absorption of Na and K has 
been shown to strongly affect the near-infrared spectra of brown 
dwarfs and extra-solar giant planets (Burrows et al. 2000). Eor 
Na and K, we used a Voigt profile V(cr-cro) in the impact region, 
up to a detuning frequency (6cr) of 30(7’/500)°® cm“' for Na 
and 50(r/500)° ® cm“^ for K, following Burrows et al. (2000). 
The Lorentz half-widths, calculated from the impact theory, are 
0.27(7’/296)"°™ cm-' atm-' forNa and 0.53(r/296)-'’ ™ cm-' 
atm"' for K. Beyond the detuning frequency, we used a profile 
in the form 


F(cr - ctq) = V{6cr){6crl{cr - (To)]^^^ 

exp[-(/ic(o- - o-o)lkT)(cr - cro)/crf], (1) 

where cr - ctq is distance from line center, V(6cr) is the Voigt 
profile at the detuning frequency 6cr, and cTf a parameter that 
we adjusted to best reproduce the absorption cross sections cal¬ 
culated by Burrows & Volobuyev (2003) for the red wings of 
the Na/K H- H 2 systems as shown in their Eig. 6 . We derived 
CTf = 5000 cm ' for Na and 1600 cm ' for K from best fitting 
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Fig. 3. Spectra for a range of models without (red) and with (blue) clouds at the distance of 10 pc and with a radius of 1 and with a radius of 

1 fJjup. For 700, 1000, 1400 and 1700 K the maximum of the blackbody is respectively around 4.1, 2.9, 2.1, and 1.7 pm 


■■■■ T.„=700K log(g) = 3.0 - T,„=1000K log(g)=4.0 . T.,^1000K log(g)=5.0 

. T,^^=1000K log(g)=3.0 - T.,,^=1400K log(g)=4.0 . T,;^1400K log(g)=5.0 

. T.„=1400K log(g)=3.0 - T,„=1700K log(g)=4.0 . T.,^1700K log(g)=5.0 

. T,^^=1700K log(g)=3.0 T.,,/=700K log(g)=5.0 - BT-SettI 1400 K 4.0 


- T.„=700K log(g)=4.0 



Temperature [K] 


Fig. 4. Temperature profiles calculated for the set of Teff and log(g) pa¬ 
rameters used to generate the spectra in Fig. 3. Triangles and squares 
on the log(g) = 4 profiles indicate the bottom of the iron and silicate 
clouds respectively. One BT-Settl profile (green) is also shown for com¬ 
parison, the star corresponding to the bottom of the dust location. The 
radiative-convective boundary is shown as a cross on the same profiles. 


of the 0.6-0.9 pm and 0.8-1.0 pm regions for Na and K, respec¬ 
tively. Profiles were calculated up to 9000 cm“^ of line center. 
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Besides line opacity, we added the collision-induced absorp¬ 
tion from H 2 -H 2 and H 2 -He using data files and subroutines pro¬ 
vided by A. Borysow^. These are based on publications by Bo- 
rysow et al. (2001) and Borysow (2002) for H 2 -H 2 , and Borysow 
et al. (1988, 1989) and Borysow & Frommhold (1989) for H 2 - 
He. 

We finally added absorption by cloud particles, discarding 
scattering. We considered condensates from Si and Fe, the two 
most abundant condensing elements in exoplanets with Teff in 
the range 500-2000 K (Lunine et al. 1989). For silicates particles, 
we used the optical constants of crystalline forsterite Mg 2 Si 04 
published by Jager et al. (2003) and for Fe liquid particles those 
from Ordal et al. (1988). 

2.2. Atmospheric model 
2.2.1. Gas composition 

The vertical profiles of H 2 O, CO, CH4, NH3, TiO, VO, Na, and 
K are calculated at each iteration from thermochemical equilib¬ 
rium, assuming a 0.83/0.17 H 2 /FIe volume mixing ratio and solar 
system elemental abundances from Table 3 of Lodders (2010). In 
the model, it is also possible to use enrichment factors over the 

^ http://WWW.astro.ku.dk/~aborysow/programs/ 
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Fig. 5. Mixing ratio profiles of important molecules assuming a cloudy atmosphere, log(g) = 4, and three different effective temperatures 700, 
1000, and 1700 K. 


solar values, independently for C, O, N, and heavier elements. 
We considered only the species that significantly affect the pro¬ 
files of the above mentioned molecules in cool giant exoplan¬ 
ets according to Burrows & Sharp (1999) and Lodders & Feg- 
ley (2006). These species are given in Table 2 (Column 2). We 
also included species involved in the formation of silicate and 
iron clouds to determine their condensation levels. Equilibrium 
abundances are derived from the equations of conservation for 
each element and using the standard molar free energies AG°(7’) 
listed in Chase (1998) to calculate equilibrium constants involv¬ 
ing the species in Table 2 (Column 2). Calculation is done level 
by level, starting from the deepest level of our grid, at highest 
pressure and temperature, and moving upward in the grid. When 
a condensate appears in a given layer, its constituent elements 
are partly removed from the gas phase and the new elemental 
abundances in the gas phase are used to calculate equilibrium 
abundances in the overlying layer. If no condensation occurs, the 
same elemental abundances are used in the overlying layer. As in 
Lodders & Fegley (2006), we take the dissolution of VO in per- 
ovskite (CaTiOa) into account, assuming an ideal solid solution 
and Henry’s law. 

2.2.2. Cloud model 

Absorption by silicate and iron clouds is included above their re¬ 
spective condensation level up to one hundredth of this pres¬ 
sure level. A particle-to-gas scale height ratio of 1 is assumed. 
We assumed spherical particles and used the Mie theory to cal¬ 
culate the absorption Mie efficiency as a function of wave¬ 
length. The particle size distribution follows a gamma distribu¬ 
tion with a mean radius r and an effective variance of 0.05. 

As discussed by Ackerman & Marley (2001) and Marley 
et al. (2012), the cloud opacity is expected to be proportional 
to the pressure pc at the condensation level, proportional to the 
total concentration of the condensing element (Si or Fe) embed¬ 
ded in various molecules at level p^, and inversely proportional 
to the gravity g. This relates to the available column density of 
condensing material at the condensation level. We write the op¬ 
tical depth of the cloud as 


Fref 

where pief = 1 bar. Because the solar elemental ratios Si/H 
and Fe/H are about the same, we assumed the particle column 
densities of the silicate and iron clouds are in the ratios of the 
pressure of their condensation levels, and thus that their Tref at 
any wavelength are in the ratios of their Qabs at this wavelength. 


We then just keep one free parameter in this cloud model, which 
is Tref for the Fe cloud at some reference wavelength. 

In running Exo-REM, we found that in some cases, when 
Tref is large and the condensation curve is close to the solution 
temperature profile, the model is instable through the iteration 
process and does not converge toward a radiative equilibrium 
solution. This is because adding particulate opacity increases 
the temperature significantly just above the condensation level. 
The temperature in this region may then become larger than the 
condensation temperature and no self-consistent solution can be 
found. This instability was also seen by Morley et al. (2012) who 
advocated a patchy atmosphere to solve this problem and reach 
a radiative equilibrium state. 

2.3. Input and output parameters 

The input parameters of the model are the effective temperature 
Teff, the acceferation of gravity g at 1 bar, which affects the atmo¬ 
spheric scafe height and thus the opticaf depth profifes, and the 
oversofar enrichment factors a for C/H, N/H, O/H and heavier 
efements X/H (a - 1 for sofar system vafues). The other set of 
free parameters are the opticaf depth of the iron cfoud Tref at 1.2 
pm and a reference condensation fevef of 1.0 bar, and the mean 
radius r of the cfoud particfes. 

For output, the modef provides the radiative-convective equi- 
fibrium temperature profife T(p), the corresponding verticaf pro- 
fifes of the absorbers at chemicaf equifibrium, and the spectrum 
at the resofution of the k-correfated coefficient distribution, i.e., 
20 cm"'. 

2.4. Examples of model outputs 

This section shows exampfes of modef outputs (spectra, temper¬ 
ature and abundance profifes) for various input parameters, af- 
fowing us to investigate the effect of surface gravity, effective 
temperature, and cfouds. Aff modefs here, unfess specified, as¬ 
sume a sofar metafficity. For modefs with cfouds (sificates and 
iron), we used Tref =1 and a mean particfe size of 30 pm (Ack¬ 
erman & Marfey 2001). We do not consider water vapor (H 2 O) 
condensation and thus formation of ice clouds, which would oc¬ 
cur in the upper atmospheres of planets with T^ff less than ~ 600 
K. 

Figure 1 allows us to compare a typical case of atmospheric 
models, with and without clouds, with two models in the liter¬ 
ature: BT-Settl and Drift-PHOENIX. The BT-Settl spectrum is 
relatively close to our case with clouds, except for some differ¬ 
ence in the wings of the alkali lines and for the presence of an ab¬ 
sorption band near 4.3 micron (outside of the wavelength range 
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accessible to SPHERE), presumably due to CO 2 , which was not 
included in our model. Other minor differences may originate 
from differences in the line lists or missing trace compounds in 
our model. The corresponding temperature profile of this BT- 
Settl model is shown in Fig. 4. On the other hand, the Drift- 
PHOENIX spectrum is very different from both the BT-Settl one 
and ours. The molecular absorption features are much less visi¬ 
ble because the clouds are thicker than in the other models and 
the spectrum ressembles that of a blackbody. 

Figure 2 presents spectra calculated for a cloud-free atmo¬ 
sphere, log(^) = 4 and two values of Teff, 700, and 2000 K. Wa¬ 
ter vapor absorbs all over the spectral range for both effective 
temperatures, with strongest bands centered at 0.94, 1.14, 1.38, 
1.87, and 2.7 pm. Other absorption features are due to other 
compounds as indicated in the figure. Besides H 2 O, K and Na 
absorption are visible both in low- and high-temperature spectra 
through their resonant lines at 767 and 589 nm, respectively. On 
the other hand, TiO, VO, and CO absorptions are only impor¬ 
tant for large Teff, while CH 4 bands play a significant role at low 
temperatures. NH 3 has a weaker effect, is only visible in planets 
with a low Teff, and provides additional absorption around 3.0, 
2.0, and 1.5 pm. 

Figure 3 shows spectra calculated for Teff varying from 700 
to 1700 K and log(g) varying from 3 to 5. The fluxes corre¬ 
spond to a planet having a Jupiter radius and located at 10 pc. 
For each set of parameters, a cloud-free model and a model with 
clouds are shown. For the case with log(g) = 5 and Tetj = 700 
K, only the cloud-free case is shown because cloud condensa¬ 
tion occurs below our pressure grid and cannot be taken into 
account. The temperature profiles corresponding to these sets 
of parameters and to a cloudy atmosphere are shown in Fig. 4. 
The locations of the iron and silicate clouds and of the radiative- 
convective boundary are also indicated. The BT-Settl tempera¬ 
ture profile for Teff = 1400 K and log(.g) = 4 (solid green line) 
is shown for comparison with the Exo-REM profile (solid black 
line). Above the radiative-convective boundary, the two profiles 
are different, with a given temperature reached one or two pres¬ 
sure scale heights deeper in the BT-Settl model. The reason for 
this discrepancy is unknown but probably lies in different verti¬ 
cal distributions of the opacity. 

For a given adding cloud absorption yields a smoother 
spectrum, decreasing the contrast between absorption bands and 
spectral windows. This is because cloud opacity, concentrated 
near the cloud base in the 1600-2100 K range, (depending on 
log(.g) and Teff, as shown in Fig. 4) more strongly reduces the 
flux in the windows, which originate from deeper levels than 
the flux in the absorption bands. Cloud opacity also affects the 
relative fluxes in the various photometric bands. Essentially, for 
a given Tefr, the flux is reduced below ~1.7 pm and increased 
longward. Therefore, the flux is significantly lower in the Y and 
J bands, at which the atmosphere is the most transparent, and 
higher in the K, L, and M bands where atmospheric opacity is 
larger. In the set of examples we show in Fig. 3, the strongest 
cloud effects are seen for log(.g) = 5 and Teff = 1000, 1400, or 
1700 K. In these cases, the emission resembles that of a black- 
body at temperature Tetj. The pressure levels where T - Teff, 
representative of the mean emission level, are the deepest of the 
set in Fig. 4, ~ 0.4 bar, and thus the cloud opacity at this level is 
the largest, according to Eq. (2). 

We only show here the effect of clouds for a single set of pa¬ 
rameters; Tief = 1 and a particle scale height equal to the pressure 
scale height. For this parametrization, the cloud optical depths at 
a given pressure level are the same for any location of the con¬ 
densation level (see Eq. 2) and any value of Tetr or log(g). There- 



Fig. 6 . Effect of varying the metallicity Z in a cloudy model with T.<r = 
1400 K and log(g) = 4.0 

fore, we cannot draw any conclusion on the relative effect of 
cloud opacity as a function of Tefr or log(g) from the calculations 
shown in Fig. 3. For example, it could be reasonable to assume 
that Tjef varies as 1/g, as does the pressure scale height, follow¬ 
ing the parametrization of Ackerman & Marley (2001, e.g., their 
Eq. 18). This would reduce the increasing effect of clouds with 
increasing gravity. Also, if the cloud is more confined near the 
condensation level than assumed here, i.e., a particle-to-gas scale 
height ratio lower than 1 , the effect of clouds would be much re¬ 
duced for cases with T^g = 700 or even 1000 K since particles 
would be confined to levels well below the mean emission level. 



Wavelength (pm) 


Fig. 7. Series of spectra calculated for a cloud-free atmosphere, T^g 
varying from 700 to 1700 K, and log(g) from 3 to 5. 

Obviously, the main effect of increasing the effective tem¬ 
perature is to increase the emitted flux but, in addition, changes 
in the spectral shape between 1 and 4 pm can be noted. As Teff 
increases, CH 4 bands, mostly visible at 1.7, 2.3, and 3.3 pm, be¬ 
come less and less intense, while the CO band at 4.7 pm as well 
as TiO and VO absorption below 1.3 pm become visible. These 
large variations may easily be detected from narrowband pho¬ 
tometry or low-resolution spectroscopy: TiO and VO signatures 
occur in the J-band, CO affects the M-band at high T^g while, at 
low Tgg, CH 4 has a strong effect in the H, K, and L-bands, and 
NH 3 has a marginal effect in the K-band. 

The spectral variations with effective temperature are, of 
course, due to changes in composition as illustrated in Fig. 5. 
Carbon is partitioned between CO and CH 4 , with a CO/CH 4 ra¬ 
tio depending on temperature and, to a lesser extent, on pressure. 
In the 700-K planet, methane dominates over carbon monoxide 
above the 1-bar level whereas in the 1700-K planet, CO dom¬ 
inates over the whole pressure grid. Similarly, nitrogen is par- 
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Fig. 8 . Maps of reduced forPictoris b SED. Vertical blue lines indicate radii of 2.5, 2.0, and 1.5 IJjup. Horizontal blue lines indicate masses of 
25, 10, 5, and 1 Mjup. From top left to bottom right, models correspond to: no cloud, t = 0.1 and < r > = 30 fim, t = 1 and < r > = 30 /rm, t = 3 
and < r > = 30 /im, t = 1 and < r > = 3 jim, respectively. 


Table 3. Clouds parameters used in the five test grids 


< r > 
(/tm) 

Aef 

'^Mg2Si04 

(A=l.2 /um; p = 1 bar) 

Tpe 

{A=\.2 /im; p = 1 bar) 


0 

0 

0 

30 

0.1 

0.015 

0.1 

30 

1 

0.15 

1 

30 

3 

0.45 

3 

3 

1 

0.018 

1 


titioned between N 2 and NH 3 , the latter being abundant in the 
observable atmosphere only for relatively low effective temper¬ 
atures (< 800 K). Also, as the effective temperature decreases, 
TiO and VO get confined to deeper levels and have thus less in¬ 
fluence on the outgoing flux. The depletion of TiO and VO in the 
upper (colder) atmosphere is due to perovskite (CaTi 03 ) forma¬ 
tion and VO condensation, respectively. Alkali Na and K affect 
all spectra in our grid but are confined at deeper levels in the 
case of low Teff atmospheres. They are removed from the upper 
atmosphere through the formation of Na 2 S condensate and KCl 
condensation, respectively. 

Fig. 6 shows the effect of varying the metallicity for given 
Teff and log(g) assuming no clouds. As expected, increasing the 
metallicity increases the depth of all absorption bands. For ex¬ 
ample, considering the water vapor band at 2.7 jum, a metallicity 
of Z = h- 0.5 produces a band depth (2.2 / 2.7 yum) twice as large 
as in the case with Z = -0.5. In principle, the metallicity of an 
observed exoplanet could thus be deduced from low-resolution 
spectroscopy provided that the temperature profile modeled from 
radiative-convective equilibrium is reliable, which also requires 
that Teff and log(g) can be accurately derived from the spectra. 
The unknown effect of clouds may be a stronger limitation in 
some cases since cloud absorption reduces the band depths and 
may mimic some decrease in metallicity. 

The gas scale height is inversely proportional to the acceler¬ 
ation of gravity g. As a result, a given optical depth at a given 
wavelength is found at deeper pressure levels when g increases. 


This explains the general behavior of the temperature profiles 
as a function of log(g) for a given effective temperature as seen 
in Fig. 4. As log(g) increases, the temperature profile generally 
moves downward along with the cloud condensation levels. The 
situation is however more complicated because of the presence 
of clouds and the dependence of molecular absorptivity with 
pressure. The effect of gravity on the calculated spectral shape is 
more subtle than that of effective temperature. It is best seen in 
spectra of Fig. 7 having no cloud opacity. Because thermochem¬ 
ical equilibrium at a given temperature depends on pressure, the 
gas abundances at a given temperature level depend on the pres¬ 
sure at this level and thus indirectly on the gravity. For example, 
the CH 4 /CO ratio at a given temperature varies as the square of 
pressure so that the methane mixing ratio at and above the atmo¬ 
spheric level where T - Teff, representative of the mean emission 
level, is larger for larger g. This explains the large increase in the 
depth of the methane bands for the Teff = 1700 K (and to a lesser 
extent 1400 K) profiles when log(g) increases from 3 to 5. In this 
case, the CH 4 mixing ratio is two orders of magnitude larger at 
the T - 1700 K level for log(g) = 5 than for log(g) = 3. These 
calculations suggest that, among objects with Teff ~ 1600-1800 
K, methane absorption would be detectable at 2.3 or 3.3 fim in 
brown dwarfs, but would probably not be detectable at 2.3 or 3.3 
/im in Jupiter-mass planets. On the other hand, for the = 700 
K profiles, the CH 4 /H 2 mixing ratio is similar for all log(g) at and 
above the T - T^g level, at its maximum value, which is twice 
the C/H elemental ratio. In conclusion, the effect of gravity on 
the spectra is significant but may be difficult to disentangle from 
compositional variations. 

3. Application of the model to actual observations 

3 . 1 . Method 

We now describe how we exploit existing data to derive charac¬ 
teristics of planets. As a first step, the model generates a grid of 
spectra for a range of physical parameters, log(g) between 2.1 
and 5.5 with a step of 0.1 and Teff between 700 and 2000 K with 
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Fig. 9. Maps of reduced for p Pictoris b J-band spectrum. From top left to bottom right, models correspond to: no cloud, t = 0.1 and < r > = 
30 pra, T = 1 and < r > = 30 pra, t = 3 and < r > = 30 pra, t = 1 and < r > = 3 pm, respectively. 


a step of 100 K. Importantly, the explored parameter space must 
be large enough to encompass all acceptable solutions. For sim¬ 
plicity, we fixed the planet radius to one Jupiter radius (Rjup), and 
leave the determination of the planet radius to the minimization 
part (see below). 

Besides the direct geometrical effect on the observed flux, 
the radius also affects the variation of the acceleration of gravity 
with altitude, and thus the scale height at a given pressure level. 
We tested this effect in a few test cases in our grid by solving 
for radiative equilibrium for two different radii and comparing 
the corresponding spectra. We only observed very small modifi¬ 
cations of the shape of the spectrum, negligible compared with 
other error bars. Hence the radius may be considered as an in¬ 
dependent scaling parameter, only affecting the observed flux 
through the area ttR^ seen from Earth. Physical parameters 
and g are derived with associated 1- or 2-cr error bars (68% and 
95% confidence level, respectively) from a analysis with n- 
1 degrees of freedom (Bevington & Robinson 2003), where n 
is the number of independent observation points (one degree of 
freedom is removed by the determination of R, see below). 

Five types of clouds were considered with characteristics 
given in Table 3: one without cloud, three with a mean parti¬ 
cle radius of 30 pm and Tref = 0.1, 1, and 3, and one with a mean 
particle radius of 3 pm and Tret = 1. The differences in the ab¬ 
sorption efficiency Qahs between iron and forsterite explain the 
differences in the optical depth calculated for each compound. 
As already mentioned in section 2.2.2, the model may be unsta¬ 
ble for Tref > 3. 

The data consist of a series of either photometric points 
(broadbands and/or narrowbands) expressed in magnitudes, a 
normalized spectrum, or both. We compute the between the 
data Jfobserved each synthetic spectrum in our grid, once in¬ 
tegrated over the photometric filters or convolved to the spectro¬ 
graph resolution (JfModei), with the following relation: 




= Z< 


Jfobserved J^Model ^2 

^2fob.served 


( 3 ) 


where AJfobserved are the uncertainties in the planet photom¬ 
etry. Then, in the case of photometric measurements, we derive 
the radius that minimizes the x^ metric 


5 logio(J^) 


i:( 


■^Observed -^Model > 


2(7 


(4) 


Therefore, the radius R (given in Rjup unit) is considered a 
global scaling parameter that does not influence the shape of the 
synthetic spectra. Finally, additional constraints based on models 
and measurements can be introduced in the analysis. Consider¬ 
ing the core-accretion model (Mordasini et al. 2012) and the hot- 
start model (Spiegel & Burrows 2012), assuming a given age of 
the star, the radius range can be restrained within lower and up¬ 
per boundaries. In addition, radial velocity measurements, when 
available, can be used to put constraints on mass and thus on g 
thanks to the relation, 


GM 


(5) 


3.2. p Pictoris b 

Table 4. Photometric measurements of Planet /3 Pictoris b. 


Filter 

Apparent Magnitude 

References 

Ys 

15.53 ± 0.34 

Males et al. (2014) 

J 

14.0 ± 0.3 

Bonnefoy et al. (2013) 

CHf s,i% 

13.18 ±0.15 

Males et al. (2014) 

H 

13.5 ±0.2 

Bonnefoy et al. (2013) 

Ks 

12.6 ±0.1 

Bonnefoy et al. (2011) 

L' 

11.02 ±0.2 

Bonnefoy et al. (2011, 2013) 

NB 4.05 

11.20 ±0.23 

Quanz et al. (2010) 

M' 

11.0 ±0.3 

Bonnefoy et al. (2013) 
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Wavelength (pm) 


Fig. 10. Top; best for P Pictoris b SED (hlack dots) with (red tri¬ 
angles) and without (hlue stars) clouds after radius and mass selection. 
Thick lines indicate the widths of the filters. CH 4 s.i% and NB 4.05 are 
very narrow filters, less than the width of plotted observation points. 
Bottom: difference between synthetic and observed fluxes divided by 
the uncertainty on the observed flux. The two horizontal lines indicate 
the ± 1 standard deviation level. 



Fig. 11. The normalized GPI J-band spechum of p Pictoris b (black) is 
plotted against the best fit of the SED also plotted in Fig. 10 with clouds 
(red triangles). Also plotted in green is a model with Teft = 1500 K, 
log(g) = 4.1 and r^f = 3, which correctly reproduces the GPI H-band 
spectrum (Fig. 12). The dotted lines represent the calculated spectra 
prior to convolution. 


For illustration, we now apply Exo-REM to the case of the 
planet p Pictoris b located at 19.44+0.05 pc (van Leeuwen 2007). 
Discovered back in 2008 (Lagrange et al. 2009), this object is a 
case study. As it orbits relatively close to its young (21+4 Myr 
Binks & Jeffries 2014) parent star, a precise followup enables a 
careful determination of the semimajor axis, which is 8.9 ^q g AU 
(Bonnefoy et al. 2014; Lagrange et al. 2013, 2014). The planet 
resides inside the circumstellar disk detected in 1987 (Smith & 
Ten'ile 1987). Importantly, Lagrange et al. (2012a) demonstrated 
that its orbital plane is in fact aligned with the warp observed by 
Mouillet et al. (1997) instead of the main disk plane, providing 
an unambiguous evidence for the disk/planet interaction. A pho¬ 
tometric event reported in 1981 could have been produced by the 
transit of this planet in front of the star (Lecavelier Des Etangs 
et al. 1997). 



800 1000 1200 1400 1600 1800 2000 


Teff [K] 



Fig. 12. Top: map of reduced for p Pictoris b H-band spectrum with 
T = 3 and < r > = 30 pm. Bottom: p Pictoris b GPI H-band spectrum 
(black dots) is compared to models with clouds with = 3. Best fit, 
in a least-squares sense, is shown as gray triangles (reff=1400 K, log(g) 
= 4.1). Models in which the gravity is changed by ± 0.3 dex are plotted 
for comparison. A model with a larger effective temperature (1500 K) 
and log(g) = 4.1, which fits the spectrum at the 1.2-cr confidence level, 
is also shown. The dotted lines represent the calculated spectra prior to 
convolution. 


We considered the whole set of available photometric mea¬ 
surements covering the near-IR wavelengths, all the way to 
the mid-IR. Observations (Table 4) were collected with NaCo 
(Lenzen et al. 2003; Rousset et al. 2003) at the VLT in the J, H, 
Ks, U, NB_4.05, M’ bands (Bonnefoy et al. 2013, 2011; Quanz 
et al. 2010) and with MagAO (Close et al. 2012) in the Ys and 
CH 4 S bands as well (Males et al. 2014). Recently, J-band (be¬ 
tween 1.12 and 1.35 pm for a resolution of 35-39) and H-band 
(between 1.51 and 1.79 pm for a resolution of 44-49) spectra 
were obtained during the GPI (Macintosh et al. 2014) commis¬ 
sioning (Bonnefoy et al. 2014; Chilcote et al. 2015). 

First, the grids of models were generated as explained here- 
above without any constraint on radius and mass. We started with 
the analysis of the photometric data points alone. The models 
with no cloud and thin clouds (Fig. 8 , top left and right panels) 
do not allow us to achieve a decent minimization, the regions of 
minima being located at the boundaries of the parameter space. 
If thicker clouds are introduced (Fig. 8 bottom), the model is 
able to reproduce the data points (reduced is lower) and the 
region limited by the 1 -cr contour falls within the grid bound¬ 
aries. We can constrain the effective temperature to 1500-1700K, 
while only a lower limit is derived for gravity (log(g) > 4). Cal¬ 
culations with < r > - 3 pm do not yield a minimum value as 
low as in the case of < r > = 30 pm and do not provide accept¬ 
able solutions at the l-cr uncertainty level (Fig. 8 bottom right). 

The same work was carried out using the J-band spectrum. In 
that case, the reduced x^ values are much lower across the grid 
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Table 5. Derived T^ff and log(g) of [} Pictoris b in each step of this analysis and shown by other studies. 


Reference 

Data type 

Constraints 

Error size 

Aff [K] 

log(g) 

Radius [Rj„p] 

Best-fit model 

This work 

SED 

no 

1 (T 

1600 ± 100 

4.0 ± 0.8 

1.68 ±0.22 

clouds, 30 fim 

This work 

SED 

no 

ICT 

1550 ± 250 

>3.2 

1.82 ±0.44 

clouds 

This work 

H-band 

no 

2cr 

1400 ± 100 

4.1 ± 0.3 

2.2 ± 0.4 

clouds, Tref = 3 

This work 

SED 

radius and mass 

1 (T 

1550 ±50 

3.8 ±0.6 

1.73 ± 0.12 

clouds, 30 fim 

This work 

SED 

radius and mass 

2o- 

1550 ± 150 

3.5 ± 1.0 

1.76 ± 0.24 

clouds 

Bonnefoy et al. (2013) 

SED 

no 

best^^ 

1700 ± 100 

4.0 ± 0.5 

1.22-1.76 

clouds. Drift Phoenix 

Currie et al. (2013) 

SED 

no 

1 cr 

1575-1650 

3.8 ±0.2 

1.65 ± 0.06 

clouds, Modified AMES-Dusty 


as a result of larger flux uncertainties (more models can fit the 
data) (Fig. 9). For models with clouds, the best correspond 
to an effective temperature of 1500 + 100 K but large gravities 
(log(g) > 3.5). In the t - 3 grid of models, the minimum 
reduced is as low as 0.14, which suggests that the error bars 
are overestimated or, more likely, strongly correlated in wave¬ 
length. With these large error bars, cloud-free models are also 
capable of fitting the spectrum at the 1 -cr level with Teff = 2000 
K, log(g) = 2.3, corresponding to a minimum reducedof 0.63. 

Assuming an age for the /? Pictoris system of 15-25 Myr, 
evolutionary models predict a radius between 0.6 and 2 Rjup- In¬ 
cluding this constraint in the x^ minimization implies a lower 
limit for the effective temperature of about 1400 K. 

Radial velocity measurements presented in Lagrange et al. 
(2012b) yield constraints on the planet mass. For separations of 
8 , 9, 10, 11, and 12 AU, the detection limit corresponds to a mass 
of 10, 12,15.5, 20, 25 Mjup, while the model-dependent mass de¬ 
rived from photometry compared to evolutionary models is > 6 
Mjup (Bonnefoy et al. 2013). Recently, Bonnefoy et al. (2014) 
used an up-to-date compilation of radial velocity measurements 
and constrained the mass limit to 20 Mjup (at 96% confidence 
level). For the purpose of being conservative, we retain a maxi¬ 
mum mass of 25 Mjup. These new constraints remove the models 
with highest gravities in the;^f^ maps. 

The mass and radius constraints allow us to more accurately 
determine the physical parameters and log(g), as presented 
in Table 5. Therefore, we propose a new determination of Teff = 
1550 + 150 K and log(g) = 3.5 + 1.0 at the 2-cr confidence level. 

These values are in good agreement with the former analyses 
by Bonnefoy et al. (2013) and Currie et al. (2013), who used 
the PHOENIX models like BT-Settl, Drift-PHOENIX or Ames- 
Dusty. 

The AMES-Cond and AMES-Dusty models represent ex¬ 
treme cases. Both solve for thermochemical equilibrium as¬ 
suming level-by-level element conservation. In the AMES-Cond 
model, all the dust is removed for opacity calculation while in 
the AMES-Dusty model the amount of dust is that derived from 
thermochemical equilibrium with no depletion process. In BT- 
Settl, the dust particle properties (number density and mean ra¬ 
dius) are derived from a comparison of timescales of various 
microphysical and transport processes. Einally, Drift-PHOENIX 
really solves for the formation and evolution of cloud particles, 
taking dust microphysics and atmospheric convection into ac¬ 
count. Bonnefoy et al. (2014) provided a detailed comparison 
with these and other models as concerns the J-band spectrum. 

Therefore, Exo-REM yields similar results as other models, 
either using photometric data or low-resolution spectra, although 
it is far less complex. A more recent study by Males et al. (2014), 
using bolometric luminosity compared to evolutionary models, 
found an effective temperature of Teff = 1643 + 32 K. In our 
work, error bars are larger probably because we consider 2 -cr 


confidence level and possibly because we probe a larger param¬ 
eter space. The models with the best;^f^ are shown in Eigs. 10 
and 11. We confirm that we need clouds to reproduce the con¬ 
trast between all photometric points. In the case with clouds we 
observe two bad fitting locations, in H band and NB4.05. As 
concerns the spectrum we have difficulty fitting the observations 
after 1.33 pm. 

Recently, Chilcote et al. (2015) published a spectrum of p 
Pictoris b in the H band (Pig. 12). We find that only models with 
thick clouds (xref = 3) can fit the data at a 3-cr confidence level. 
Note that published error bars of the spectrum only account for 
random errors (at the 1 - 2 % level) and do not incorporate sys¬ 
tematic uncertainties, such as an estimated 10 % uncertainty on 
the overall spectral slope. Our best fit, in a least-squares sense, 
is obtained for T^g - 1400 K, log(g) = 4.1 and = 3 (Pig. 12) 
with a reduced of 1.12. However, this case would imply a 
too large radius of 2.1 Rjup to reproduce the absolute flux, which 
is based on the photometric measurement of Males et al. (2014). 
We also show a case with T^g- 1500 K, log(g) =4.1 and Tref = 3 
yielding = 1.18, which agrees within error bars with the pa¬ 
rameters derived from photometric measurements (Table 5) with 
mass and radius constraints and from the GPI J-spectrum (Pig. 
11). The agreement of both synthetic spectra with the GPI H- 
spectrum is very good. In Table 5, we give the model parameters 
that fit the spectrum at the 2 -cr level and not for the 1 -cr level be¬ 
cause, as mentioned above, the error bars of Chilcote et al. may 
miss some systematic uncertainty. Adding a 10% uncertainty on 
the overall spectral shape would significantly increase the error 
bars on the derived parameters. The fits we obtain are closer to 
the observations than models presented in Chilcote et al. (2015). 
The differences between the models possibly result from differ¬ 
ent modeling of the cloud opacity. To fit the shape of the ob¬ 
served spectrum, models with thick clouds are needed. Models 
with no or thin clouds produce too much contrast between the 
peak and both ends of the observed spectrum. In Pig. 12, we 
also show cases where log(g) is varied by 0.3 dex around the 
previous case to illustrate the strong sensitivity of the shape of 
the spectrum to this parameter (as gravity increases the spectrum 
broadens). Spectral observations in this band thus provide a way 
to constrain the gravity of the planet, although one would need 
to investigate to what extent it can be disentangled from cloud 
opacity and metallicity. 


4. SPHERE expected observations 

In this section we discuss the ability of SPHERE to put useful 
constraints on gravity and effective temperature according to the 
quality of the data. The purpose is to link the photometric er¬ 
rors to the uncertainties in the physical parameters of planetary 
atmospheres. A related analysis to derive log(g) and T^g was per¬ 
formed by Vigan et al. (2010) with the narrowband differential 
filters of SPHERE combined with AMES-Cond/Dusty (Allard 
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et al. 2001, 2003), BT-Settl (Allard et al. 2007) and Burrows 
models (Burrows et al. 2006). 

In the following, we consider 12 test cases for which the 
model has a robust convergence, and which cover a represen¬ 
tative range of log(g) = 2.5, 3.5, 4.5, - 800, 1100, 1400, 

1700 K, and cloud properties (xief = 1 for 30 fim particles). For 
this preliminary analysis, we focus on the near- IR broadband fil¬ 
ters Y, J, H, Ks, offered in IRDIS, the SPHERE camera, as well 
as the Y-H mode of lES (39 wavelengths, with 0.014-0.020 jum 
between adjacent pixels), the near-IR spectrograph (Table 6). 

The spectra of test cases were integrated over IRDIS filters 
and a photometric error was added to the integrated flux to mimic 
an actual photometric measurement. We considered several error 
amplitudes (in magnitude). Am = 0.01, 0.05, 0.1, 0.5, and 1.0, 
corresponding to very good to very poor data. The same error 
amplitude was applied to all filters, and we did not consider data 
with various qualities. In fact, the flux error in the SPHERE mea¬ 
surements is expected to show some correlation between wave¬ 
lengths that affect the spectral shape to some extent. We did not 
take these systematic errors into account in the present analysis. 
As in the previous section, the photometry of test cases was com¬ 
pared to the grid of models using the minimization (no mass 
or radius constraint). The uncertainties in physical parameters, 
ATeff and Alog(g), were derived from the 2-cr contours. Results 
are shown in Eig. 13. We observe that ATeff decreases as log(g) 
increases and conversely, Alog(g) increases as T^g increases. 

When photometric errors are small, on the order of 0.01 mag, 
the errors are often smaller or similar to the step of the grid. 
At the other extreme, when the photometric error is as large as 
1 mag, all models contained in the grid match the observation, 
hence the errors on the physical parameters exceed the range of 
the grid. We conclude that the effective temperature and grav¬ 
ity can be constrained to 200 K and 0.5 dex, respectively, if an 
accuracy of 0.2 mag is achieved. 

We also considered flux spectra, normalized to unity at the 
peak, with an error of 0.01, 0.05, or 0.1, constant for all wave¬ 
lengths. The same exercise performed with our set of synthetic 
spectra (Eig. 14) indicates accuracies of 200 K for T^g and 0.5 
dex for g, assuming a precision of 0.1. 

The number of available photometric data points, as well as 
the covered spectral range, also have an impact on the accu¬ 
racy of the retrieved physical parameters. Eor instance, consid¬ 
ering the combination of two broadband filters, the set H-hKs is 
more appropriate to constrain the effective temperature (Eig. 15 
a) while for gravity a large spectral range is preferable (such as 
J-tKs; Eig. 15 b). We now consider three possible sets of pho¬ 
tometric data points: two SPHERE filters (H and Ks, or J and 
Ks), four SPHERE filters (Y, J, H, Ks), and finally the same four 
SPHERE filters with two NaCo filters (L’, M’). We assume an 
accuracy of 0.1 mag on SPHERE data. To achieve the same ac¬ 
curacy in the physical parameters as previously achieved (200 
K for Tgg and 0.5 dex for g), we conclude that at least three 
data points are required. In addition, a significant improvement 
is achieved if the SPHERE photometry is complemented with 
the NaCo MIR filters. With L’ and M’ filters we can expect un¬ 
certainties in Teff lower than 100 K and also a smaller error on 
log(g). 


5. Conclusions 

To analyze photometric and spectroscopic data from new 
instruments like SPHERE at the VET we have developed 
Exo-REM, a radiative-convective equilibrium model to simulate 
the atmosphere of young Jupiters, which are privileged targets 
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Fig. 13. Effect of uncertainties in the magnitude of photometric data 
points upon uncertainties in derived Teff (top) and log(g) (bottom). 
Cases with 2-cr error bars exceeding our test grid are not plotted. 


Table 6. Characteristics of SPHERE IRDIS filters and lES spectro¬ 
scopic mode 


Name 

Central wavelength [/tm] 

EWHM* [fim] 

BB Y 

1.0425 

0.139 

BB J 

1.2575 

0.197 

BBH 

1.6255 

0.291 

BB Ks 

2.1813 

0.3135 

Y-H 

0.957-1.636 

R~30 


★: Eull width at half maximum 


for direct imaging of exoplanets. The model incorporates 
opacity from the molecules and atoms that are relevant to 
observable levels for giant exoplanets having T^g < 2000 K. 
It assumes that vertical profiles of these species are governed 
by thermochemical equilibrium. Cloud absorption by iron and 
silicate clouds is included through a simplified formalism and a 
limited number of free parameters. 


Article number, page 11 of 15page.l5 















A&A proofs: manuscript no. Baudinoetal_17aout 


500 


400 


: 300 


c 

o 


200 - 


100 - 


0 

Uncertainties on normalized spectrum 


O 800K, log(g)=2.5 


■ 


700 



O 800K, log(g)=2.5 

[> HOOK, log(g)=2.5 






HOOK, log(g)=2.5 ' 

0> 1400K, log(g)=2.5 







1400K, log(g)=2.5 

1700K, iog{q)=2.5 







^ 1700K, log(g)=2.5 

^ 800K, log{g}=3.5 

o 

♦ . 


600 

♦ 


H 800K, log(g)=3.5 

^ HOOK, log(g)=3.5 





^ HOOK, log(g}=3.5 

A 1400K, iog(g)=3.5 


► 





^ 1400K, log(g)=3.5 

1700K, loa(a)=3.5 


, _ , 

500 

■ # 


1700K, log(g)=3.5 

tt 800K, log(g)=4.5 





tt 800K, log(g)=4.5 ' 

HOOK, iog(g)=4.5 



'—' 




HOOK, log(g)=4.5 

' A 1400K, log(g)=4.5 






^ 1400K, log{g}=4.5 

1700K, log(g)=4.5 




400 

. ► 

♦ 

1700K, log{g)=4.5 ^ 



C 









o 





■ 

► 


L- 

o 

L. 

300 

■ > 


■ 




LU 

200 




, 








i 




100 

0 


► 

- m-m- 

1_■_ 

o.bi 


0.1 


2 

4 

6 


Number of photometric points 



Uncertainties on normalized spectrum 



Fig. 14. Effect of uncertainties in the normalized spectra upon uncer¬ 
tainties in derived Teif (top) and log(g) (bottom). Cases with 2-cr error 
bars exceeding our test grid are not plotted. 


Fig. 15. Effect of number of points in the SED upon uncertainties in 
derived I’eff (top) and log(g) (bottom). The two-point case corresponds 
to the J and Ks filters. Cases with 2-cr error bars exceeding our test grid 
are not plotted. 


We used Exo-REM to analyze data available for Pictoris b 
and to derive physical parameters of the planet. We inferred an 
effective temperature Teff = 1550 ± 150 K, log(g) = 3.5 + 1, and 
a radius R - 1.76 + 0.24 Rj^p (2-cr error bars) from photometric 
measurements and considering independent constraints on mass 
and radius. These results are similar to those previously derived 
by other authors using different atmospheric models. The differ¬ 
ence is that we considered 2-cr error bars (rather than 1-cr) and 
explored a wider range of parameters, in particular, with lower 
values of g, than in previous studies. Our 2-cr uncertainties in¬ 
clude measurement error as well as model dependence on our 
limited set of cloud parameters (optical depth, particle radius). 

We were also able to reproduce the H-spectrum of Chilcote 
et al. (2015) within their (small) error bars in contrast to other 
models displayed in that paper. Using this spectrum alone, the 
derived parameters are; Tetj = 1400 + 100 K, log(g) = 4.1 + 0.3, 
and R = 1.9 + 0.2 Rjup (2-cr error bars). 

We investigated the ability of SPHERE to characterize 
exoplanets with the IRDIS broadband filters and the Y-H 
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spectroscopic mode of lES. A couple of filters (H, Ks) appear 
best suited to constrain Tes, while the couple (J, Ks) is more 
appropriate to constrain log(g). Combining MIR NaCo L’ and 
M’ observations with SPHERE photometry enables us to obtain 
good constraints on both T^ff and log(g). 

We plan to explore the set of free parameters of Exo-REM 
more systematically than shown in this paper. In particular, 
we will more extensively study the effect of metallicity and of 
cloud parameters (scale height, reference optical depth, particle 
size). In future works, we may consider constraints coming from 
ab-initio models like BT-Settl or Drift-Phoenix, which provide 
guidelines for a range of realistic and physical cloud parame¬ 
ters. We will also add absorption by water ice particles that are 
expected to form in giant exoplanets having lower effective tem¬ 
peratures than studied here. 

SPHERE was commissioned successfully and the instrument 
is open to the community. Known planets are prime targets for a 
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thorough characterization (Lagrange et al. in prep.; Zurlo et al in 
prep.). 
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Appendix A: Radiative-convective equiiibrium 
modei, numericai method 

In a one-dimensional radiative-convective equilibrium model, 
the net flux (radiative -H convective) is assumed to be constant 
as a function of pressure level. This net flux ttF is equal to 

nF = (A.l) 

where is the effective temperature of the planet. We first 
solve for purely radiative equilibrium and neglect heating from 
the parent star. This is justified as long as we are interested in 
hot young giant exoplanets relatively far from their parent star. 
Assuming a planet with - 700 K at a distance of 7 AU from 
the star, the stellar flux absorbed by the planet would amount to 
less than 0.1% of the planet’s thermal emission. 

Discarding scattering, the net flux at pressure level p, in a 
plane parallel geometry, is given by 


OO OO 


0 tAp) 


7tF(p) do- J" B^{t'^)E 2 {t'^ - Ta-)dT'^- 

tAp) 

f 


BAT'o-)E2{T<r - T^dr'^ 


, (A.2) 


where To-(p) is the optical depth at pressure level p and 
wavenumber cr, is the Planck function at the temperature 

of level of optical depth and wavenumber cr, and £2 is the 
second-order exponential integral. 


The integral over wavenumber is calculated over the range 
[c^min, CTniax], with cTmin = 20 cm^' and cTniax = 16000 cm^^ 
This range is sliced into Wo- intervals of width So- = 20 cm"\ 
over which the Planck function is taken as constant. The radia¬ 
tive transfer integral over each interval of width So- is calculated 
through a correlated-k distribution method with quadrature 
points (Goody & Yung 1989). The atmospheric grid consists of 
Hp atmospheric levels equally spaced in ln(p) between pressure 
levels pnjax at the bottom of the grid (J = 1) and rip = 64). As¬ 
suming a linear variation of the Planck function B with optical 
depth T within any layer [pj-i, pj], i.e. 

BAr^) = BAt^j-i) + BAT^j) ~ , (A.3) 

"^crj -1 crj t'crj -1 T (j-.j 

the contribution of this layer to the flux at wavenumber cr can be 
analytically calculated and expressed as a linear combination of 
the Planck functions at pressure levels pj and pj+i. We also add a 
contribution from below the atmospheric grid (p > Pmax) assum¬ 
ing a semi-infinite layer with the same variation of the Planck 
function as in Eq. ( A.3 ) for the first layer. Summing over all 
layers and spectral intervals, the net flux at level pj can then be 
expressed as 

n,r "p np 

TrFipj) - 2nSo- zz ^<r.(7j')^wiA(cri,pj,pj,,0, (A.4) 

i=i j'=i 1=1 

where Ba-XTy) is the Planck function at temperature of the /* 
pressure level of the grid and wavenumber cr; at the middle of 
the spectral interval of width So-. The parameter vix is the 
weight applied to the Z* set of k-correlated coefficients used in 
the quadrature for the spectral integration over any spectral inter- 

np 

val (Z wi = !)■ The parameter A(cri,pj,pj/, 1) is a dimensionless 
1=1 
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factor that couples pressure levels pj and py and only depends 
on the grid of optical depths for the Z* set of k- correlated coef¬ 
ficients of the Z* spectral interval. 

We then search for the temperature profile that ensures radia¬ 
tive equilibrium, i.e., 

mpj) - o-Tt^ (A.5) 

for j varying from 2 to Wp. We do not use the flux at the first, 
deepest level as a constraint because the variation of the Planck 
function at deeper levels is fixed arbitrarily in the model to that 
of the first layer. To solve this system of np_i equations, we use 
a constrained linear inversion method described in Vinatier et al. 
(2007) and based on Conrath et al. (1998). The algorithm mini¬ 
mizes the quadratic difference (x^) between desired (o-T^A and 
calculated fluxes with the additional constraint that the solution 
temperature profile lies close to the reference profile. Starting 
from an initial guess profile Tq, an approximate solution Ti is 
derived from the equation 

T„ = + aSK^C-^AF, (A.6) 

with n - I, where AF is the difference vector between the de¬ 
sired and calculated fluxes (crTefr^ - nF(pj)), K is the kernel 
matrix with Ajj- equal to the derivative of the flux at level pj 
with respect to the temperature at level py, 5 is a normalizecl 
two-point Gaussian correlation matrix that provides a vertical 
filtering of the solution needed to avoid numerical instabilities, 
and a a scalar parameter that controls the emphasis placed on 
the proximity of the solution Ti to the reference profile Tq. We 
used a correlation length of 0.4 pressure scale height. The kernel 
matrix is calculated from Eq. (A.4), neglecting the dependence 
of A with temperature, which is generally much weaker than that 
of the Planck function, i.e., 

^ dBAATy) ^ 

Ajj- = InSo- ^^ zrriA(cri, Pj, Pj,, Z). {A.l) 

i=i j' 1=1 

Matrix C is equal to 

C^aKSK^ + E, (A.8) 

where £ is a diagonal matrix with £jj/ equal to the square of the 
flux error acceptable at the f' pressure level, usually set to 0.1% 
of o-Tl^. 

The nonlinearity of the problem requires an iterative process 
in which £„ is obtained from Eq. (A.6) after updating the ref¬ 
erence profile to Tn-i and recalculating the kernel matrix K for 
profile Tn-i. The iteration process is pursued untilis less than 
1 and no longer significantly decreases. The a parameter in Eq. 
{A.l) is chosen to be small enough to ensure convergence and 
large enough to reduce the number of iterations needed. Typi¬ 
cally ten iterations are needed. Note that the final solution does 
not depend on the initial profile Tq or on the choice of a. Eor 
To, we used one of the three temperature profiles calculated by 
Allard et al. (2003) for = 900, 1300, and 1700 K. We choose 
that having Teft closest to the input value to ensure rapid conver¬ 
gence. 

In a second step, the solution profile is checked against 
convective instability by comparing the model lapse rate V = 

fni/ip^'i) adiabatic value = RjCp, where R is the gas 

constant and Cp the temperature-dependent molar heat capac¬ 
ity for the H 2 -He atmosphere. Regions where the lapse rate ex¬ 
ceeds the adiabatic value are unstable against convection. They 
are found in the bottom of the pressure grid p > pad, with pad 
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being the level where the lapse rate equals the adiabatic value. In 
that case, convective heat transfer occurs setting back the lapse 
rate to the adiabatic value. In Exo-REM, we do not solve ex¬ 
plicitly for convection. We add a convective flux to the radiative 
flux in Eq. (A.4) through an analytical function that is essentially 
zero when V < Vad and rapidly gets very large when V > Vad. 
We found that the following function: 

Eeonv(Pj) = lO-Vr^V^OOtV/V,.-!] (^ 9 ) 

is adequate to ensure negligible superadiabaticity in the final so¬ 
lution profile. We then set the lapse rate of the purely radiative 
solution to the adiabatic value, plus a small amount (0.015), at 
levels p > Pad ■ The kernel matrix K is calculated adding 

gEcon^v(Pj) ^ Vad fa(Pj/Pj-l)) (a) 

to Aijj and 

dF (p,) ^ Vadfa(Ej/;^j-i)) ib), 

Ol\-\ 

(A. 10) 

to /qj_i in Eq.(A.7). The iterative process is finally restarted with 
the modified flux and kernel until convergence is achieved. Typ¬ 
ically another set of ten iterations is needed. Our model in this 
paper has 64 pressure levels equally spaced in ln(p) between 50 
bar and 0.01 mbar. 
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